#include <stdlib.h>
#include <stdio.h>
#include "nbody.h"
#include <math.h>
#include <string.h>
#include <pthread.h>
#include <time.h>
#define MAX 255

int numCuerpos, numPasos, anchuraTablero, alturaTablero, numProcesos;
double g, intervalo, minDist;
double *posicion[2];
double *velocidad[2];
double *fuerza[2];
double *masa;
int *indices[4];
pthread_mutex_t mp;

int main(int argc, char **argv){
    int i=0;
    Hilos *hilo;
    
    numProcesos=(atoi(argv[3]));
    leerCabecera(argv[1]);
    calcularIndices();
    //  pthread_mutex_init(&mp,NULL);
    for(i=0;i<numPasos;i++){
        crearHilosCalcFuerzas(hilo);
        crearHilosCalcPosicion(hilo);
    }
    printf("Guardando datos en el fichero: %s\n", argv[2]);
    escribirDatos(argv[2]);
    
    return 1;
}

int leerCabecera(char *fichero_entrada){
    FILE *file;
    int  cont;
    
    //Apertura del fichero y lectura de la cabecera
    
    file = fopen(fichero_entrada, "r");
    
    fscanf(file, "%d %lf %lf %d %lf %d %d", &numCuerpos, &g, &intervalo, &numPasos, &minDist, &anchuraTablero, &alturaTablero);
    
    inicializarDatos();
    
    printf("numcuerpos: %d, g: %lf, intervalo: %lf, numPasos: %d,mindist: %lf,anchura: %d,altura: %d\n", numCuerpos, g, intervalo, numPasos, minDist, anchuraTablero, alturaTablero);
    
    for(cont=0; cont<numCuerpos; cont++){
        fscanf(file, "%lf %lf %lf %lf %lf \n", &posicion[0][cont], &posicion[1][cont], &velocidad[0][cont], &velocidad[1][cont], &masa[cont]);
        
        
    }
    return 1;
}
    
    int escribirDatos(char *archivo){
        FILE *file;
        int cont;
        
        //Apertura del fichero y escritura de la cabecera
        file = fopen(archivo, "wr");
        fprintf(file, "%d %lf %lf %d %lf %d %d\n", numCuerpos, g, intervalo, numPasos, minDist, anchuraTablero, alturaTablero);
        
        //Escritura de los datos
        for(cont=0; cont<numCuerpos; cont++){
            fprintf(file, "%lf %lf %lf %lf %lf\n", posicion[0][cont], posicion[1][cont], velocidad[0][cont], velocidad[1][cont], masa[cont]);
        }
        fclose(file);
        return 1;
    }
    
    void inicializarDatos(){
        posicion[0] = (double *)malloc(numCuerpos* sizeof (double));
        posicion[1] = (double *)malloc(numCuerpos* sizeof (double));
        fuerza[0] = (double *)malloc(numCuerpos* sizeof (double));
        fuerza[1] = (double *)malloc(numCuerpos* sizeof (double));
        velocidad[0] = (double *)malloc(numCuerpos* sizeof (double));
        velocidad[1] = (double *)malloc(numCuerpos* sizeof (double));
        masa = (double *)malloc(numCuerpos* sizeof (double));
        indices[0]=(int *)malloc(numProcesos*sizeof(int));
        indices[1]=(int *)malloc(numProcesos*sizeof(int));
        indices[2]=(int *)malloc(numProcesos*sizeof(int));
        indices[3]=(int *)malloc(numProcesos*sizeof(int));
    }
    
    void *calcularFuerzas(int proceso){
        double fax, fay, dist, f, alpha, aux;
        int i, j, planetaInicial, cuerposProceso;
        int fin=0;
        
        cuerposProceso=(numCuerpos/numProcesos);
        planetaInicial=cuerposProceso*proceso;
        
        
        for (i = 0; i < numCuerpos&&proceso==0; i++) {
            fuerza[0][i] = fuerza[1][i] = 0;
        }
        
        for (i = indices[0][proceso]; (i <= indices[2][proceso]&& !fin); i++) {
            
            if (i== indices[0][proceso])
                j=  indices[1][proceso];
            else
                j=  i+1;
            
            for (; j<numCuerpos && !fin; j++) {
                
                dist = sqrt((posicion[0][j] - posicion[0][i]) * (posicion[0][j] - posicion[0][i]) + (posicion[1][j] - posicion[1][i]) * (posicion[1][j] - posicion[1][i]));
                
                if (dist > minDist ) {
                    
                    f = g * masa[i] * masa[j] / (dist * dist);
                    
                    if (f > MAX_F)
                        f = MAX_F;
                    
                    aux  = (posicion[1][j]-posicion[1][i]) / (posicion[0][j]-posicion[0][i]);
                    
                    alpha = atan(aux);
                    
                    fax   = f * mycos(alpha);
                    fay   = f * mysin(alpha);
                    
                    pthread_mutex_lock(&mp);
                    fuerza[0][i] += fax;
                    fuerza[1][i] += fay;
                    fuerza[0][j] -= fax;
                    fuerza[1][j] -= fay;
                    pthread_mutex_unlock(&mp);
                    
                }
                if((j+1)==indices[3][proceso]&&i == indices[2][proceso])
                    fin=1;
            }
        }
        pthread_exit((void*)0);
    }
    
    void *calcularPosicion(int proceso){
        double *ax, *ay;
        int i, planetaInicial, cuerposProceso;
        
        cuerposProceso=(numCuerpos/numProcesos);
        planetaInicial=cuerposProceso*proceso;
        
        ax = (double *)malloc(sizeof(double) * numCuerpos);
        ay = (double *)malloc(sizeof(double) * numCuerpos);
        
        for (i = planetaInicial; i < planetaInicial+cuerposProceso; i++) {
            
            //calculo de las aceleraciones
            ax[i] = fuerza[0][i] / masa[i];
            ay[i] = fuerza[1][i] / masa[i];
            //calculo de las velocidades
            velocidad[0][i] = velocidad[0][i] + ax[i] * intervalo;
            velocidad[1][i] = velocidad[1][i] + ay[i] * intervalo;
            //calculo de las posiciones
            posicion[0][i] = posicion[0][i] + velocidad[0][i] * intervalo;
            posicion[1][i] = posicion[1][i] + velocidad[1][i] * intervalo;
            // printf("[%d] px: %lf vx: %lf intervalo: %lf\n",i,posicion[0][i],velocidad[0][i],intervalo);
            if (posicion[0][i] <= 0) {
                posicion[0][i] = 2;
                velocidad[0][i] = - velocidad[0][i];
            }
            if (posicion[0][i] >= anchuraTablero ) {
                posicion[0][i] = anchuraTablero - 2;
                velocidad[0][i] = - velocidad[0][i];
            }
            if (posicion[1][i] <= 0) {
                posicion[1][i] = 2;
                velocidad[1][i] = - velocidad[1][i];
            }
            if (posicion[1][i] >= alturaTablero) {
                posicion[1][i] = alturaTablero - 2;
                velocidad[1][i] = - velocidad[1][i];
            }
            //printf("body %d: x=%d, y=%d vx=%lf vy=%lf ax=%lf ay=%lf fx=%lf  fy=%lf\n", i, x[i], y[i],vx[i],vy[i],ax[i],ay[i],fx[i],fy[i]);
        }
        free(ax);
        free(ay);
        pthread_exit((void*)0);
    }
    
    void calcularIndices(){
        
        int numIteracionesTotales=0, numIteProceso, i, j, aux=0, procesoActual=1;
        
        for(i=1; i<numCuerpos; i++)
            numIteracionesTotales+=i;
        
        numIteProceso=numIteracionesTotales/numProcesos;
        
        indices[0][0]=0;
        indices[1][0]=1;
        
        for (i=0;i<numCuerpos&&procesoActual<(numProcesos);i++) {
            for(j=i+1;j<numCuerpos&&procesoActual<(numProcesos);j++) {
                aux++;
                if(aux==numIteProceso) {
                    
                    if((j+1)<numCuerpos) {
                        indices[2][procesoActual-1]=i;
                        indices[3][procesoActual-1]=j+1;
                        indices[0][procesoActual]=i;
                        indices[1][procesoActual++]=j+1;
                    }
                    else {
                        indices[2][procesoActual-1]=i+1;
                        indices[3][procesoActual-1]=i+2;
                        indices[0][procesoActual]=i+1;
                        indices[1][procesoActual++]=i+2;
                    }
                    aux=0;
                }
            }
        }
        
        if((j)<numCuerpos) {
            indices[0][numProcesos-1]=i-1;
            indices[1][numProcesos-1]=j;
        }
        else {
            indices[0][numProcesos-1]=i-1;
            indices[1][numProcesos-1]=i+1;
        }
        indices[2][numProcesos-1]=numCuerpos;
        indices[3][numProcesos-1]=numCuerpos;
        
    }
    
    int crearHilosCalcFuerzas(Hilos *hilo) {
        int i=0;
        hilo=(Hilos *)malloc(numProcesos*sizeof(Hilos));
        for (;i<numProcesos;i++) {
            hilo[i].idHilo=i;
            
            if(pthread_create(&hilo[i].thread, NULL, calcularFuerzas,  hilo[i].idHilo))
                printf("Error en la creacion del hilo [%d]\n", hilo[i].idHilo);
        }
        for (i=0;i<numProcesos;i++)
            pthread_join(hilo[i].thread, NULL);    
  
        return 1;
    }
    
    int crearHilosCalcPosicion(Hilos *hilo) {
        int i=0;
        hilo=(Hilos *)malloc(numProcesos*sizeof(Hilos));
        for (;i<numProcesos;i++) {
            hilo[i].idHilo=i;
            
            if(pthread_create(&hilo[i].thread, NULL, calcularPosicion,  hilo[i].idHilo))
                  printf("Error en la creacion del hilo [%d]\n", hilo[i].idHilo);
     }
        for (i=0;i<numProcesos;i++) 
            pthread_join(hilo[i].thread, NULL);
        return 1;
    }
    double factorial(double x) {
        double fac = 1;
        int i;
        for (i=2; i<=x; i++) fac *= i;
        return fac;
    }
    
    double mypow(double x, double exp) {
        double result = 1;
        int i;
        
        result = 1;
        for (i=0; i<exp; i++) result *= x;
        return result;
    }
    
    double mysin(double alpha) {
        double result = 0;
        int i;
        
        for (i=0; i < MAX_ACCURACY;i++) {
            result += (double)(mypow(-1, i) / factorial((double)(2*i+1))) * (double)mypow(alpha, 2*i+1);
        }
        return result;
    }
    
    double mycos(double alpha) {
        double result = 0;
        int i;
        
        for (i=0; i < MAX_ACCURACY;i++) {
            result += (double)(mypow(-1, i) / factorial((double)(2*i))) * (double)mypow(alpha, 2*i);
        }
        return result;
    }

